Non-equilibrium phase transitions in biomolecular signal transduction 
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We study a mechanism for reliable switching in biomolecular signal-transduction cascades. Steady 
bistable states are created by system-size cooperative effects in populations of proteins, in spite of 
the fact that the phosphorylation-state transitions of any molecule, by means of which the switch is 
implemented, are highly stochastic. The emergence of switching is a nonequilibrium phase transition 
in an energetically driven, dissipative system described by a master equation. We use operator and 
functional integral methods from reaction-diffusion theory to solve for the phase structure, noise 
spectrum, and escape trajectories and first-passage times of a class of minimal models of switches, 
showing how all critical properties for switch behavior can be computed within a unified framework. 



I. INTRODUCTION 

A. The emergence of devices from biomolecular 

systems 

Two related questions define a fundamental role for 
statistical physics in systems biology: (1). How do 
biomolecular systems achieve reliable "device-level" be- 
havior when they consist of highly stochastic compo- 
nentry [l|, in particular molecular complexes held to- 
gether by low-energy hydrogen or Van der Waals bonds? 
(2). Are such systems structured in a modular fash- 
ion i.e. can complex molecular networks decom- 
pose into quasi-autonomous functional units performing 
identifiable tasks which are robust against many physical 
parameter changes and are recombinable in evolution? 

The device logic of biomolecular systems employs 
many of the same abstractions as electronic engineer- 
ing [3j, including amplification, filtering, and switching. 
Switches (which use elements of amplification and filter- 
ing) are used wherever a discrete sequence of events is 
required, enabling committed threshold response to en- 
vironmental cues during development [U, H[ , establishing 
"checkpoints" for intermediate states, in processes like 
cell division that require strict sequencing of events @, 
or programming the complex progressions of cell-type dif- 
ferentiation [2]. 

Two classes of biomolecular processes in which modu- 
lar switching is widely recognized are gene expression and 
signal transduction. In gene expression, the switch states 
are associated with patterns of genes activated for protein 
production, and a stochastic component is introduced in 
the system by the small number of weakly-bound tran- 
scription factors [7] . In signal transduction, the states of 
the switch are frequently determined by the number of 
phosphate or methyl groups covalently attached to spe- 
cific target amino acid residues on one or more dedicated 
proteins §. Stochasticity in this system can again be 
caused by the small number of proteins or the complexes 
these proteins form with the catalysts which promote 



attachment or detachment of the phosphate or methyl 
groups. Gene expression changes cell type on slower 
timescales (minutes to hours) than signal transduction 
(seconds), but since many cell- type changes occur in re- 
sponse to external signals [1, [B[ , or rely on amplification 
and stabilization of internal signals Q , switching behav- 
ior is often produced by both systems acting together. 



In this paper we consider the problems of mechanism 
for the emergence of robust switching in stochastic molec- 
ular systems, and of quantitative estimation of the noise 
and stability properties of switches. We consider switch- 
ing in signal transduction via the mechanism of phospho- 
rylation (addition of phosphate groups via the action of 
a kinase) and dephosphorylation (removal of phosphate 
groups by the action of a phosphatase) since these are 
a common motif found in most if not all signalling net- 
works. In addition, the relative simplicity of phosphory- 
lation transitions lends itself to the abstraction of many 
real transduction cascades in terms of a few processes, 
allowing us to isolate the problems of formation, con- 
trol, and robustness of the switch. We observe that the 
most fundamental unit in signal-transduction cascades is 
a single type of target protein with multiple phospho- 
rylation states (called phosphoepitopes), among which 
transitions are naturally modeled as a reaction-diffusion 
process. The network properties necessary for switching, 
which may be distributed in real systems among several 
proteins in a signal-transduction cascade, or between the 
cascade and the genetic transcription factors for the tar- 
get proteins, are readily lumped together and assigned to 
a single species to produce minimal models. This coarse- 
graining reduces network complexity while still keeping 
its essential regulatory features. We are able to write 
down exact master equations for such ideal models, and 
to solve them systematically with field-theoretic methods 
from reaction-diffusion theory. 
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B. Senses of "switching", their uses, and how they 
are achieved robustly 

Switching in biomolecular systems, at the least, refers 
to sigmoidal response to input signals, termed "ultrasen- 
sitivity" Q . Such response implies a sharp sigmoidal but 
continuous response in the concentration of a molecule 
over a narrow range of a (stationary) signal. In con- 
trast, bistability pj| is a form of switching made pos- 
sible when two stable states, Si and S2, co-exist over 
a signal range. As a consequence, bistable systems ex- 
hibit two distinct thresholds as the signal is varied, one 
at which a transition occurs from Si to S2 and another 
at which the system switches back from S2 to Si. The 
separation of thresholds leads to path dependence or hys- 
teresis, and makes a switched state more impervious to 
stochastic fluctuations of the signal around the transi- 
tion point, than in the ultrasensitive case. Hysteresis 
may be obtaine d by adding positive feedback to sig- 
moidal response Weak hysteresis may be used to 
stabilize input signals, equivalent to signal "debouncing" 
in engineering, while strong hysteresis leads to bistabil- 
ity, toggling, and long-term memory [l2|, fl3| - We will 
look specifically at toggling, because this subsumes all 
the other phenomena of interest. 

The molecular mechanisms used repeatedly in signal 
transduction for amplification, threshold sensitivity, and 
switching, are shown in Fig. [TJ as they are instantiated 
in the Mitogen- Activated Protein Kinase (MAPK) family 
of transduction cascades. This widely duplicated and di- 
versified homologue family is used througout the eukary- 
ote kingdom [14j . mostly to regulate gene expression in 
response to cell-membrane received signals. MAPK cas- 
cades employ three proteins, each with an unphosphory- 
lated state, and respectively one, two, and two phospho- 
rylated states. Phosphorylation and dephosphorylation 
of each protein is catalyzed by exogenous kinases and 
phosphatases, and in addition the fully-phosphorylated 
states of the first two proteins act as kinases (phosphory- 
lation catalysts) on the proteins following them in the 
cascade. The cascade is thus an actively driven, dis- 
sipative system, maintained away from equilibrium by 
the supply of activated (high-energy) phosphate donors. 
When used for switching, MAPK cascades may have pos- 
itive feedback from the output to the highest-level pro- 
tein, either through gene expression or through inhibition 
of degradation of the active state @ . 

The structure of MAPK and other cascades was ab- 
stracted by Goldbeter and Koshland H to the minimal 
system shown in Fig. [51 which they propose as the sig- 
naling counterpart to the transistor (a better analogy 
would be to the bistable flip-flop, as they use it). A single 
protein species has a single phosphorylation site. Phos- 
phorylation and dephosphorylation occur via the action 
of catalysts/enzymes with which the protein can form 
enzyme-substrate complexes. Depending on the rate of 
these reactions, the steady state fractions of the phospho- 
rylated (or unphosphorylated) protein can vary abruptly 




FIG. 1: The catalytic topology of the mitogen-activated pro- 
tein kinase (MAPK) cascades (from Ref. with dashed 
positive feedback arrow representing polyadenylation that in- 
hibits degradation of the top protein). Each level in the dia- 
gram represents the phosphorylation states of a single protein, 
with phosphorylation and dephosphorylation transitions rep- 
resented by arrows. Exogenous and internal catalysis is rep- 
resented by (other) arrows pointing to the transition arrows. 



as a function of these rates. This analogy has been ex- 
tended to an elaborate analysis of the properties [T3| and 
combinatorial logic Q of such switches. In particular, 
12] considers the case where the the phosphorylated epi- 
tope acts as an intermolecular autocatalyst on phospho- 
rylation transitions of any unphosphorylated proteins in 
the population, and shows that this leads to bistability. 
However, autocatalytic feedback only creates a bistable 
switch if the response of the underlying phosphorylation 
chain is sigmoidal fill ], which in this model requires sat- 
uration of the exogenous phosphatase rate via the forma- 
tion of catalyst-substrate complexes as an intermediate 
step between the unphosphorylated and phosphorylated 
states of the protein. 
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FIG. 2: The Goldbeter-Koshland minimal phosphoryla- 
tion couple. A single protein species has an unphosphory- 
lated state W and a singly-phosphorylated state W* . Ref. [8| 
chained such couples to model intermolecular phosphoryla- 
tion, in combination with exogenous catalysts that are often 
present for both transitions as well. 



We note, however, that kinetic control through satu- 
rated complex formation is not the only way to obtain 
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sigmoidal response, because with two or more phospho- 
rylation sites per protein, the concentration of the fully- 
phosphorylated state is a sigmoidal function of the ra- 
tio of exogenous kinase to phosphatase, even when all 
catalysts act in the linear proportional regime (in other 
words, when catalyst-substrate complexes act to catalyze 
transitions effectively instantaneously, and are limited 
only by their frequency of formation through binary en- 
counters). This occurs as long as the catalytic activity 
is distributive i.e., if at most one modification (phos- 
phorylation or dephosphorylation) takes place at each 
enzyme-substrate encounter [l5j . and ordered (if succes- 
sive phosphorylations take place at different residues, an 
ordered mechanism implies that dephosphorylation takes 
place in strictly the inverse order) [To ]. 

Two of the MAPK proteins have this structure, and 
more significantly, the intermolecular catalysis within 
the cascade is nonspecific to phosphorylation reactions 
on a given protein, though each transition is catalyzed 
through an independent event @. We show below 
that, combining this form of sigmoidality with positive 
feedback, it is possible to obtain bistability through a 
non- equilibrium phase transition, in which the individ- 
ual events of catalysis leads to a polarized distribution 
of phosphorylation states of the target protein. Such 
population-level cooperative effects, (proposed also in the 
context of genetic switches in Q), bestow the stability of 
macroscopic (thermodynamic) systems on the otherwise 
highly stochastic events of phosphorylation and dephos- 
phorylation. We suggest that the properties of phase- 
transition-mcdiatcd switching are one source of adaptive 
preference for multiple phosphorylation sites and non- 
specific catalysis, which one encounters repeatedly (his- 
tidine kinase cascades may have as many as 26 phospho- 
rylation sites (l7j)- 

Previous studies have also shown that multisite phos- 
phorylation with saturation kinetics at each modification 
step can lead to bistability even in the absence of feed- 
back [II GJ]. Hence both kinetic control and population- 
level polarization can lead robustly to bistability in some 
parameter domains. However, the two mechanisms are 
distinguished by their responses to mutations and by 
their control parameters. Single-molecule control causes 
switching properties to change if rate kinetics change, in 
a way that population polarization does not, while the 
role of nonspecific catalysis in models with population- 
level cooperative effects (at least in the form we will con- 
sider) creates a different kind of sensitivity. An important 
constraint on the evolutionary innovation, preservation, 
and diversification of phenotype (any expressed function- 
ality) is the shape of its neutral network [2(J Hl[ (the 
degenerate space in the genotypc/phcnotype map with 
regard to that functionality). The phenotype of phase- 
transition-mcdiatcd switching is more nearly controlled 
by the topology of the catalytic network than by its ki- 
netics, an idea that has been proposed as a source of 
robustness in the segment-polarity network Q, and the- 
oretically grounded in the case of general enzyme-driven 



reaction networks in [l9| . 

Note that we do not study spatio-temporal correla- 
tions induced by diffusion of enzymes and / or enzyme in- 
activation. A recent study 22[ shows that even with 
the enzymes acting according to a distributive mecha- 
nism, rapid rebindings of the enzyme molecules to the 
substrate molecules can lead to a loss of ultrasenstivity 
and bistability. We do not consider the effect of protein 
degradation either. Our model is however a first theoret- 
ical fully stochastic study of the MAPK cascade, mod- 
elled earlier, to our knowledge, only via rate equations 
[1 E H HI Hi or stochastic simulations H, EH]. 

In this context, we study quantitatively the three criti- 
cal properties of a phase-transition mediated switch: the 
conditions for existence of bistability, the noise charac- 
teristics of those fluctuations that preserve the domain 
in the bistable phase, and the large excursions that limit 
memory or reliability of the switch, and which near the 
threshold for bistability, can lead to finite-particle num- 
ber corrections to that threshold. These have only been 
considered piecemeal before in other models, with condi- 
tions for bistability treated in the infinite-particle (deter- 
ministic) limit [ToL H3] , noise from internal and external 
sources related through ad hoc response functions [28$ . 
and stability treated at the level of bounds on scaling, 
for systems already assumed reduced to one relevant di- 
mension [13j. 



C. Reducing to appropriate models 

Most biological literature on this subject focuses on 
phenomenological modeling of (usually mean-field behav- 
ior in) observed or designed systems @, 0-H i, 0, E3 ■ 
We are however more interested in the possibility of sta- 
tistically motivated universality classification of strate- 
gies for switching, which might explain evolutionary reg- 
ularities in cascades. Therefore, in addition to idealiz- 
ing molecular mechanisms responsible for sigmoidal re- 
sponse and positive feedback as properties of single pro- 
tein species, to make the polarization-based equivalent of 
the flip-flop from adding feedback to our Fig. [2] (equiva- 
lent to Fig. 12 of Ref. Q), we advisedly exploit symme- 
try, of either the catalytic topology or the parameters, 
to make analysis tractable. This approach also aids in 
decomposing effects responsible for switching, and relat- 
ing these to other equilibrium or non-equilibrium phase 
transitions. Thus our minimal models deliberately differ 
from the familiar cascade families in areas not directly 
related to the production of switching [2^]. The model 
of Markevich et al [l8[ is also in this category in demon- 
strating bistability (via kinetic control) at the level of a 
single stage of the MAPK cascade. 

The model we propose for a cooperative- 
phosphorylation switch is shown in Fig. [3j Each of 
N molecules of a single type of protein has J phospho- 
rylation sites indexed j G 1, . . . , J, which we suppose for 
simplicity to be phosphorylated and dephosphorylated 
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in a definite order. All phosphorylations are catalyzed 
by exogenous kinases, and all dephosphorylations by 
exogenous phosphatases. Because the enzymes are 
assumed to operate in the linear regime where complex 
formation is not rate limiting, the catalytic rate per 
reaction is proportional to the numbers I and P of 
kinase and phosphatase particles respectively (We set 
the constant of proportionality equal to one by choice of 
the units of time). The site modifications occur in a spe- 
cific order, thus sidestepping combinatorial complexity 
Furthermore, phosphorylation and dephosphorylation of 
substrate proteins is assumed to follow a distributive 
mechanism, whereby a kinase (phosphatase) enzyme 
dissociates from its substrate between subsequent mod- 
ification events [13, 31]. Hence the substrate has J + 1 
states. 



a b 




FIG. 3: Multisite phosphorylation and feedback structure, 
(a) This panel depicts the basic phosphorylation chain with- 
out feedback in which a target protein with J sites is phospho- 
rylated by a kinase / and dephosphorylated by a phosphatase 
P. The ordered succession of phosphorylations yields J + 1 
modification states, labelled 0, 1, . . ., J. (b) The fully phos- 
phorylated target protein relays a signal into pathways that 
eventually feed back on the phosphorylation chain, (c) Simpli- 
fication of (b) in which the fully phosphorylated target protein 
acquires kinase activity and directly feeds back on the chain. 
We refer to this network configuration as the asymmetric cir- 
cuit, (d) Schematic of the network with symmetric feedback 
in which the substrate protein is bifunctional, whereby the 
fully (de)phosphorylated form catalyzes (de)phosphorylation 
of its own precursors. 

We obtain positive feedback from intermolecular au- 
tocatalysis; specifically, proteins in the j — J state are 
kinases that act interchangeably with the exogenous ki- 
nases (unequal catalytic power between j = J and exoge- 
nous kinases can easily be added, at the cost of another 



parameter). The phosphorylation chain with feedback is 
shown in the bottom half of Fig. [3] Panel (c) depicts 
an asymmetric topology in which the fully phosphory- 
lated substrate catalyzes its own phosphorylation, while 
panel (d) shows the symmetric version in which a sub- 
strate molecule is bifunctional, acting as both a kinase 
and phosphatase depending on its modification state. Ki- 
nase / and phosphatase P are exogenous forces on the 
modification of the substrate, but the feedback is an en- 
dogenous force whose strength is proportional to the oc- 
cupancy of the end-states of the chain. This occupancy is 
subject to intrinsic fluctuations and depends on the total 
number of substrate molecules. 

The assumption of intermolecular autocatalysis is stan- 
dard [13] , and we consider below the self-consistent back- 
grounds with kinase-only autocatalysis, as the nearest 
equivalent to the kinetically-controlled switch Q. In or- 
der to analyse the model beyond mean-field theory (using 
field theoretic techniques) we go further in the interest 
of simplicity, and symmetrize the topology as in Fig. [3ji, 
by making the unphosphorylated state (indexed j = 0) a 
phosphatase, interchangeable with the exogenous phos- 
phatases. 

The feedback topology of the model caricatures a few 
elements present in biological systems. One such ele- 
ment is the competition between antagonistic pathways 
that may underlie cellular decision processes (for exam- 
ple [Hj]). A multisite phosphorylation chain of the type 
considered here could function as an evaluation point be- 
tween competing and antagonistic pathways influenced 
by different active phosphoforms of the chain, provided 
these pathways feed back to the chain. In a less extreme 
case, the fully phosphorylated form activates another ki- 
nase which then interacts with the chain. In these sce- 
narios, feedback is mediated by a series of intervening 
processes, which may well affect the propagation of fluc- 
tuations. Yet, if delays are not too large, the collapsed 
scheme of Fig. |3t could be a reasonable proxy with the 
added benefit of mathematical tractability. 

A scenario corresponding more literally to our model 
involves a bifunctional substrate capable of both kinase 
and phosphatase activity, depending on the substrate's 
modification state. One example is the HPr kinase/P- 
Ser-HPr phosphatase (HprK/P) protein, which operates 
in the phosphocnolpyruvatexarbohydrate phosphotrans- 
ferase system of gram-positive bacteria. Upon stimula- 
tion by fructose-l,6-bisphosphate, HprK/P catalyzes the 
phosphorylation of HPr at a seryl residue, while inorganic 
phosphate stimulates the opposing activity of dephospho- 
rylating the seryl-phosphorylated HPr (P-Ser-HPr) [33] ]. 
Another example of a bifunctional kinase/phosphatase is 
the NPJI (Nitrogen Regulator II) protein. It phospho- 
rylates and dephosphorylates NRI. NRI and NRII con- 
stitute a bacterial two-component signaling system, in 
which NRII is the "transmitter" and NRI the "receiver" 
that controls gene expression. NRII autophosphorylates 
at a histidine residue and transfers that phosphoryl group 
to NRI. The phosphatase activity of NRII is stimulated 
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by the PII signaling protein (which also inhibits the ki- 
nase activity). Several other transmitters in bacterial 
two-component systems seem to possess bifunctional ki- 
nase/phosphatase activity [34j . 

Both the network with asymmetric topology (auto- 
kinase only, Fig. (3J;) and the network with symmetric 
topology (Fig. [5J1) but asymmetric catalytic concentra- 
tions I P undergo formally first-order phase transi- 
tions, so that regions of bistability are always metastable 
at finite N. However, in the topologically symmetric 
case, these continue smoothly through a second-order 
transition at / = P, in which symmetry of both topol- 
ogy and parameters ensures exact bistability with finite 
residence time in domains, at all N where the phase tran- 
sitions exists. This simplification permits us to estimate 
the residence times with an expansion in semiclassical 
stationary points of an effective action, without encoun- 
tering the complexities of path integrals for metastable 
processes [35[ , though numerically we expect this also to 
be a good approximation to residence times in metastable 
states with similar "barrier heights" in the first-order 
case. Symmetry also permits the closed-form computa- 
tion of the noise kernel about the monostable phase with 
a unique equilibrium, which generates a natural measure 
for "weakness" or "strongness" of the first-order tran- 
sitions at nearby values of I/P as a function of N. We 
therefore perform a thorough analysis of the second-order 
transition, to establish methods and provide a reference 
solution to qualitatively understand the mechanisms of 
bistability and metastability in the more general cases 
with similar stochastic structure. 



D. Methods of treatment for the stochastic 
problem 

While differential equations for mean chemical concen- 
trations (the current standard method of analysis) can 
give good estimates of the existence of hysteresis and 
bistability when approximating systems with as few as 
tens to hundreds of molecules, they of course preclude 
the treatment of noise, fluctuation-induced corrections 
to mean-field behavior at small particle number or near 
critical points, and large excursions such as domain flips 
(when the system switches from one bistable state to an- 
other). Pure mass-action models also ignore spatial con- 
straints such as scaffolding by the cytoskeleton or the 
proteins themselves, and the dimensionality of physical 
diffusion in the cytosol or membranes. 

A better approximation is given by the master equa- 
tion for the probability of instantaneous particle distri- 
butions in models like that of Fig. [31 which in princi- 
ple captures all orders of stochastic processes, though 
such simple models still omit spatial effects. The general 
properties of the master-equation (in the diffusion, or 
"Fokker-Planck" approximation) for a one-dimensional 
switch have been used to obtain scaling relations and 
loose bounds on the stability achievable from such a 



switch as a function of the number of molecules it em- 
ploys 0. 

Operator methods, analogous to Hamiltonian methods 
in quantum mechanics, have been developed in reaction- 
diffusion theory (see 36] for a review) to efficiently 
handle the collective excitations that diagonalize general 
master equations without time-reversal symmetry. These 
have been used in the context of gene expression [3] to 
estimate the number of stable cell types made possible by 
many randomly combined transcription factors, making 
use of similarities to ground states of random-bond Ising 
models. 

From the operator-valued evolution kernel, one can 
obtain an equivalent path-integral representatio n by ex- 
panding at each time in a basis of coherent states (37l.[38|. 
Stationary-field expansion in the path integral gener- 
alizes the classical differential equation for concentra- 
tions to consistently incorporate fluctuation effects (by 
means of a perturbatively-corrcctcd effective action 39] ) , 
and the sum over "approximate stationary points" of lo- 
cally least-action identify the typical configuration his- 
tories associated with domain flips. More sophisti- 
cated approaches, similar to those used here, have also 
been used to incorporate fluctuation effects into efficient 
lumped-parameter expansions for networks with multiple 
timescales (ioj . 

Master equations can also be solved numerically by the 
Gillespie algorithm [ill ]. or simulated directly, and we 
use such simulations to validate our anlaytic results be- 
low. The lack of convenient symmetries in real biomolec- 
ular systems promises to make analysis intractable for 
most quantitative phenomenology, and recourse to nu- 
merics is likely to be the only general-purpose solution. 
However, the path integral's separation of moments in a 
natural small-parameter expansion, and of perturbative 
noise from formally non-perturbative large excursions, 
provides an intuitive decomposition of the mechanisms 
fundamental to switching and stability. At mean-field ap- 
proximation, we find surprising similarities of the phase 
transition in this driven system to the magnetization 
transition in the discrete, equilibrium, mean-field Ising 
ferromagnet, and a transition between this classical crit- 
ical behavior at finite J and a condensation effect more 
similar to Bose-Einstein condensation at J — > oo. The al- 
gebraic distinction between self-consistent backgrounds, 
perturbative fluctuations, and non-perturbative domain 
flips, elementary in the analysis, is also a subtle distinc- 
tion, difficult to make without systematic measurement 
biases, in the numerics. 



E. Main results from the analytical treatment 

The mean-field results, which are reported in detail 
in [29| . and which can also be recovered from our 
effective-action treatment in this paper, reproduce the 
standard differential equations for mass action. The 
stationary states arise from conditions of detailed bal- 
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ance between phosphorylation and dephosphorylation, 
self-consistent with the concentrations they produce of 
autocatalytic phosphoepitopes in relation to exogenous 
catalysts. Specifically, we show how both symmetric and 
asymmetric topologies create domains of mono- and bi- 
stability in the parameter space (I/N,P/N), and how 
the population asymmetry in the self-consistent state de- 
pends on the coupling g = N/y/ IP and exogenous asym- 
metry I /P. 

The perturbative expansion in Gaussian fluctuations 
about the self-consistent background provides a system- 
atic construction of the noise spectrum of the phos- 
phorylation chain. At lowest order it predicts a cusp 
~ 1/ \g — g c | in the variance of the order parameter, 
equivalent to the Curie- Weiss prediction for the spin-1/2 
mean- field ferromagnet. More surprising, we find that 
the entire perturbative approximation to the noise spec- 
trum on all sites is generated from a single bare mode, 
effectively coupled to a single Langevin field. This result 
replaces the ad hoc noise kernels one must entertain in 
the absence of a first- principles treatment [ID, H2| ■ 

The nonperturbative expansion in semiclassical con- 
figurations of locally least action predicts the leading 
large- N dependence of the domain residence time in the 
bistable regime, as a function of the dimensionless rates 
of the problem I/P and V IP /N (though here we solve 
only for the symmetric case I/P = 1, where bistabil- 
ity remains exact at finite particle numbers). These 
configurations, the dissipative equivalent to the instan- 
tons of Euclidean equilibrium field theory |35j , solve two 
problems. First, from the high-dimensional configuration 
space of the A^-particle, (J + l)-site chain, they extract 
the one-dimensional contour of most likely configurations 
to mediate domain flips, assumed given in Ref. [l3| . 
Second, the action along this trajectory, oc N at fixed 
(I/N, P/N), is the leading exponential in the residence 
time, for which Ref. [l3[ correctly predicts the scaling but 
gives no algorithm to compute the coefficient (known in 
large-deviations literature as the rate function [43[). 

Similar leading-exponential dependencies have been 
computed in Ref's. [42|,|44j. For reference to this work, 
we note that the passage to the diffusion limit or Fokker- 
Planck equation in Ref. [l3|], and the closely- related 
use of the Gaussian approximation for fluctuations in 
Ref. [42], [55] are formally uncontrolled approximations, 



whose limitations and ranges of validity are pointed out 
in Ref. [4J]. One purpose for our paper is to present the 
larger systematic analysis within which such approxima- 
tions arise. 



F. Layout of the paper 

Sec. |TT] introduces the master equation for the model 
class of Fig. c and d, and derives the phase diagram 
for steady states from conditions of detailed balance of 
the mean particle numbers. Sec. IIIII converts the master 
equation, first into the equivalent representation in terms 
of a state in a Hilbert space, and then into the equivalent 
path-integral representation through an expansion in in- 
termediate Poisson distributions. Sec. lIVI derives the per- 
turbative expansion in fluctuations about the mean fields 
of the path integral, including the equivalent representa- 
tion in terms of a Langevin equation, and the leading- 
order perturbative approximation to the fluctuations in 
the order parameter. Sec. |V]then considers the enlarged 
expansion in approximate stationary points needed to de- 
rive the trajectories and rate of domain flips. Finally, 
Sec. I VII summarizes the consequences of these technical 
results for the conceptual understanding of biomolecular 
signal transduction and switching. 



II. MASTER EQUATION AND MEAN-FIELD 
BACKGROUNDS 



An instantaneous configuration of N proteins on the 
,7+1 sites of Fig. [3] defines a vector n = (no, n\, .. . , nj), 
where rij is the number on site j. Fixed particle number 
implies that n lives on the integer lattice in the J-simplex 



N. We denote a (generally time-dependent) 



probability distribution on configurations P(n), and sup- 
press the time index t in the notation. 

A stochastic process for particle hopping is completely 
defined by the master equation for P(n), which is the 
"probability inheritance" equation induced by the tran- 
sition probabilities on the simplex. For Fig. [3] with cat- 
alytic rates proportional to the number of catalytic par- 
ticles, this is 



j-i 



dt 



+ (P + n Q - S Qj ) (n j+1 + 1) P{n - lj + - (P + n ) n j+1 P(n)] . 

I 



(1) 



where lj denotes the vector with jth component equal Our assumption that phosphorylation and dephospho- 
to 1 and all other components zero. rylation happen in a definite sequence makes transition 
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rates from site j proportional to rij and the catalyst con- 
centration, without additional combinatorial factors. For 
the asymmetric (auto-kinase only) topology, the factors 
no and Sq j in the second line of Eq. ([T]) are absent, and 
dephosphorylation depends only on the rij and the ex- 
ogenous phosphatase number P. 

Time-dependent average particle numbers on each site 
are defined as 



(2) 



and it is easy to see from Eq. {T]) that Ylj=o d ( n j) /dt = 
0. 

It is also useful to write the equation for the center-of- 
mass of the system defined as C = J2 n j 3 n jP{ n )- This 
becomes 



d_ 

di 



C = N (I - P) + [P (n a ) - I (nj)} 
+ N[(nj)- (no)] + [(r%)-(n*j)] 



(3) 



As we will see later, this exact equation can be used to 
estimate the fluctuations of the order parameter for large 
J. 

In what follows, we first look at the mean-field approx- 
imation already elaborated on in (29[. The mean-field 
approximation for evolution of (rij) under Eq. ([1} re- 
places all joint expectations with products of marginals: 
(n-juj) w (rij) (nj), etc. 



A. Detailed balance: symmetric topology 

Under the mean-field approximation, the system of N 
interacting particles essentially decouples into a system 
of N independent particles executing a random walk on 
the lattice of J + 1 sites. Within this approximation de- 
tailed balance of phosphorylation and dephosphorylation 
between adjacent sites in the chain holds, and depends 
on a catalytic ratio which we will denote x. For the 
symmetric topology, x = (I + (nj)) / (P + (no)), and we 
recognize two convenient nondimensional parameters: 



and 



N 

Zip 



a- 



(4) 



(5) 



As autocatalysis, scaled by N, induces bistable order (i.e. 
it favours configurations in which most particles are piled 
up towards one or the other end of the chain), and exoge- 
nous catalysis, scaled by I and P, induces homogeneity 
(configurations in which particles are uniformly spread 
out on the chain), g is the coupling strength of the model, 
with strong coupling favoring broken symmetry. I/P is 
then the measure of exogenous asymmetry. 



In terms of these and the fractional occupations 
(rij) /N, the catalytic ratio may be written 



+ g ( nj ) / N 

X — nr. ; , ttt = . 



(6) 



e-V 2 + g (n ) /N 
By induction on j, time-independent solutions satisfy 

(n d ) = at (n Q ) , < j < J, (7) 
and the normalization 



• 7 i _ r J+i 

N = J2(n j ) = (no)^— ■ 

From Eq. ^ and Eq. (J7J we can evaluate 

2sinh(i^)sinh(^) 

iV <nj - ?l ° > = sinh(^) 

and we can rewrite Eq. (|SJ) as 



(8) 



(9) 



\nj - np) 
N 



2sinh(|) sinh(^) 



sinh(^±ie) 



(10) 



When Eq. (j9J) and Eq. (fT0| are nonzero, they have the 
ratio 



sinh(^) sinh(^±i^) 



sinh(|) sinh(^i^) 



(11) 



For I ^ P, Eq. (fTTj) always holds (though it may be 
negative or singular), while for I = P we have the possi- 
bility of the degenerate case where (nj — no) = and g is 
unconstrained. For g < g c (\ = 0) (a second-order criti- 
cal point) this is the stable asymptotic distribution, while 
for g > gc(^ — 0) it is unstable. The graph of g versus 
£ for a few (non-positive) values of A at J = 2 is shown 
in Fig. [4] (Positive A generate curves reflected through 
£ = 0.) The graph defines a pseudo- inverse £(<?), which 
gives the stationary solutions within the mean-field ap- 
proxiation. Where £ is triple-valued (not a well-defined 
inverse), the central branch is in all cases unstable, and 
the two outer branches are stable. 

The character of the curves in Fig. 2] is preserved for 
all J > 2, though the derivative of the stable curve for 
I = P above its critical point becomes discontinuous at 
£ = for J — > oo. This discontinuity is related to the 
transition from Curie- Weiss to Bose-Einstein-like behav- 
ior of the order parameter, discussed below. The curves 
corresponding to all A have regular limits at large J. 

We can identify a set of <?c(A) as the local minima in 
Fig. @] above which the £(g) graph becomes triple- valued. 
The A — > limit of these minima smoothly converges on 
the second-order critical coupling 



g c (X = 0) 



■7+1 
J— 1' 



(12) 



8 



15 



10 



TV 

, \ \ 
\\\' 

•\\\\ 

\\\\v 

', \ \ \ \ \ 
'. \ \ \ \ 

V\V\\ N 

- \ \ \ *. \ 

- \ \ \ \ \ 

w . \ \ \ \ 
> \ \ \ s s 



\ 



FIG. 4: g(£) from Eq. along contours of constant A, for 
which the two branches proceed outward from the £ = line 
in order of increasing |A| (different A values are depicted by 
different line styles). The branches with the lowest and largest 
values of A are marked on the figure. The other A values can 
be read on the left from the intersection with the abscissa 
where f = A. Where <?(£) has a single-valued inverse, that 
function £(<?) defines the unique steady-state distributions. 
Where the inverse is triple- valued, the largest |£| are stable 
solutions, and the central branch is unstable. 



Converting the pair (A,5c(A)) to (I/N,P/N) values 
yields the phase diagram shown in Fig. [5] for a range of J 
values. The interior region I ~ P and sufficiently small 
\/ IP /N = 1/g is bistable, and outside this region the 
sign of £ = log a; equals that of A = log I /P. As we 
demonstrate in [29|], these theoretical estimates match 
very well with data from Monte carlo simulations. 



B. Detailed balance: asymmetric topology 

In the auto-kinase-only asymmetric model, positive A 
(I > P) is never bistable, because the particles are al- 
ready biased toward rij, the only site with positive feed- 
back. Therefore we graph only A < 0, though the alge- 
braic solutions are valid everywhere. 

Instead of Eq. ©, the catalytic ratio is x = 
(I + (nj)) /P, which reduces in nondimensional param- 
eters to 



_ e x ' 2 +g(nj) /N _ 



e -A/2 - • 

Equations (|7|) and © still hold, but instead of Eq. 
we choose the reduction 



(13) 



N 



2 sinh 



exp[(J-i)£] 



(14) 




X J = 99 



<i=9 



'i =5 



0.2 0.3 0.4 0.5 

l/N 



0.7 0.8 0.9 



FIG. 5: Phase diagram for the symmetric topology, from the 
minima gc(A) of Fig. [4] and their equivalents for a range of 
J. The regions inside the chevrons are bistable. 



The appropriate reduction of Eq. 
Eq. ([TO]), is now 



dSJ, counterpart to 



N 



sinh 



(I) 



exp(^) sinh(^) 



(15) 



Eq. (fl~4| and Eq. (fT5|) are regular at all ^, so we always 
have a defined function g(f), of the form 



2 sinh S inh(i±l^) 
sinh (|) exp 



(16) 



A graph at J = 2, which is the asymmetric-topology 
counterpart to Fig. SI is shown in Fig. [Bl In the bistable 
phase, there are still three branches for £ at given g, 
with the outer two stable and the central one unstable. 
The obvious differences are that now there is a maximal 
A for bistability, and that the leftmost stable branch at 
any A moves positively in £ as g increases because of the 
asymmetric topology, whereas in the symmetric topology 
it moved negatively in £. 

At any A below a (negative) J-dependent threshold, we 
can extract the minimal and maximal g values for bista- 
bility (below the minimum, £ follows A = log (I/P) qual- 
itatively; above the maximum, only the largest-^ branch 
is stable because of too-strong positive feedback). In- 
verting 5 m in(A) and (A) to (I/N,P/N), we obtain 
the phase diagram for bistability shown in Fig. [71 

The upper boundary of each bistable region, defined 
by 3min(A), is a distorted counterpart to the upper gcW 
branch in Fig. [SJ and the two converge to the same limit 
as J — > oo (where feedback from no becomes irrelevant). 
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FIG. 6: <?(£) along contours of constant A for the asymmetric 
topology, J = 2 and the same A values as in Fig. [4] The 
maximal A for bistability generates the curve whose minimum 
derivative is zero. 



formed this continuation smoothly by weighing the no 
catalytic strength with a parameter e € [0,1].) Fur- 
ther, the first-order transition in g along any I/P ray 
in the symmetric topology continues smoothly through 
the second-order transition at I/P = 1, at the apex of 
the domain of bistability. Whether the first-order transi- 
tions in the neigborhood of I/P = 1 are strong or weak 
depends on the ^-support of the stationary solution for 
P(n) (a function of N), in relation to the difference be- 
tween the stable mean values at gc (A). 



C. Phase transition and order parameter versus J 

We now restrict attention to the case of symmetric 
topology and exogenous catalysis setting I = P = q, and 
consider the behavior of the natural mean-field order pa- 
rameter \{nj — no) \ /N as a function of J. Expanding 
Eq. (fTU|) in small £, and inverting Eq. (fTTj) relative to 
g c = gc(^ = 0) in Eq. (|12l) . we find the mean- field crit- 
ical scaling of the discrete Ising ferromagnet, up to a 
J-dependent prefactor: 




J = 99 



0.8 0.9 



FIG. 7: The phase diagram for auto-kinase feedback only. 
The domains of bistability are somewhat smaller than in 
Fig. [5] and are shifted toward smaller I/P, but otherwise 
the characteristics are qualitatively and even quantitatively 
similar. 



The lower boundary, defined by g m ax(A), replaces the 
reflected lower gc(A) branch in Fig. [5j and converges to 
the diagonal / = P at J — > oo. 

Thus we see that classically, the first-order phase tran- 
sitions are similar for symmetric and asymmetric feed- 
back topology, one being deformable into the other in 
the (I/N,P/N) parameter space. (We could have per- 



\{nj - n ) | 
N 



J + i Vs. 



1/2 



(17) 



The small-£ approximation is valid for g — g c < g c — 
1, above which the order parameter saturates to a J- 
independent envelope value 



\{nj - n ) | 
N 



1 

1 - -. 

9 



(18) 



The exact mean-field prediction for \{nj — rao)| /N versus 
g from Equations (fTU)) and (TTTT) is compared to numerical 
simulations for J + 1 = [5, 10, 100], in Fig. [5] 

Since g c — 1 — > 2/ J for large J, Eq. (fP8"|) also gives the 
behavior in the formal J — > oo limit. The derivative of 
the order parameter converges to one in arbitrarily small 
neighborhoods of the critical point, rather than to oo as 
in the Curie- Weiss regime; thus J — > oo defines a dif- 
ferent universality class than any finite J. Qualitatively, 
the distinction between small and large- J is determined 
by whether one or both reflecting boundaries are sensed 
by the near-critical symmetry-broken state. The large- 
J transition resembles Bose-Einstein condensation in the 
sense that either no or nj accounts for a finite fraction of 
the particles, with the remainder "thermalized" with an 
exponential distribution in j into the interior, at "tem- 
perature" self-consistently determined by q/N — l/g. To 
understand the nature of these transitions beyond mean- 
field theory we introduce the operator and path-integral 
representations of the master equation and its solutions. 
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FIG. 8: Order parameter \(rij — no)\ /N for the symmet- 
ric phosphorylation chain: mean-field theory (lines) and sim- 
ulations (symbols). J + 1 = [5,10,100] corresponds to 
[dot, dash, solid] for lines and [+,o, x] for symbols. Parti- 
cle numbers used in the simulations are respectively N = 
[4000, 2000,400]. 



III. OPERATOR, STATE, AND PATH 
INTEGRAL REPRESENTATIONS 

A. Operators, states, and time evolution 

The operator representation of master equations from 
reaction-diffusion theory [3(| |38[ begins by introducing 
raising and lowering operators aj, aj for each site on 
the lattice, with the commutation relations of orthogonal 



a,-, a. 



6 jf . These 



quantum harmonic oscillators 
define a Hilbert space through their action on a "vacuum" 
state dj |0) = 0, Vj, and its conjugate (0| a - = 0,Vj. 

Number states indexed by the vector n are defined 
through the action of the raising operators 



A classical distribution P{n) has the state representa- 
tion in the n basis 



(22) 



\ip) is equivalent to a generating function f(z) of a 
( J + l)-component complex vector z, under the associ- 
ation O Zj, aj -f-> d/dzj. Glauber normalization is 
equivalent to a prescription for shifting Zj — > Zj + 1, and 
evaluating the resulting function at Zj ~ 0,Vj. A thor- 
ough treatment of these methods for handling generating 
functions and functionals is provided in Ref. |45| in the 
context of the analysis of master equations for evolution- 
ary games. However for the treatment that follows below, 
we need only the definitions provided above in order to 
proceed. 

The master equation (fTJ) corresponds to state evolution 
equation known as the Liouville equation 



dt 



|V) = -«IV0- 



(23) 



in which the nonlinear, diffusive Liouville operator that 
evolves the state in time is given by 



.7-1 

n = 9 ( a i+i 



H a j+ i - H 

1 J V Q 



(24) 

Here / = P = q as mentioned earlier. The differential 
equation (|23p is formally reduced to quadrature to give 
the time-dependent state relation 



-fit 



(25) 



Normalization of P(n) and the number states \n) implies 
(0| exp a^j \ip t ) — 1, Vt. We further recognize the ex- 
ogenous catalytic strength q = 1/t as defining a natural 
timescale, and the natural coupling g = N/q in Eq. 
as before. 



n (»])"• io). 



(19) 



and differ from quantum-mechanical number states in be- 
ing normalized with respect to a universal Glauber state 



(0|e 



= l.Vn. 



(20) 



B. Coherent-state expansion and path integral 

At weak nonlinearity (small g), it is both intuitive and 
computationally efficient to expand solutions to Eq. (125[) 
in eigenvectors of the annihilation operators aj [38j |. 
which are the Poisson distributions in rij . We start with a 
normalized initial state arbitrarily parametrized by mean 
occupation numbers 



a]a 3 , 



The number operator for each j is defined as fij 
and extracts the appropriate coefficient from n in the 
Glauber norm, 



| exp 




(21) 



IV'o) = exp 



J> a}-l 



10) 



(26) 



in which judicious choice of the fij cancels surface terms 
associated with transients. (Self-consistency of these pa- 
rameters with stationarity under il may be used from 
the operator representation to obtain moments of P(n), 
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as was done in Ref. Q, though we will proceed directly 
to the time-dependent field action here.) To form a basis 
for coherent-state expansion (again, see 38] for detials 
of this procedure) at increments of time, we introduce a 
complex-valued vector field 4>t = (4>0i <t>i> ■ ■ ■ , 4>j)t > an< ^ 
its adjoint <fr\. At a set of t = kAt, we insert the repre- 
sentation of identity 



h' -4>\ 



' |0) (0|e« 



E 



n) (n| 



into 

through 
(0|ex P (E : 



(0|exp 
Eq. (USD 



and 



Eq. 

- 1 



(27) 
expressed 
(J2BJ) as 

^ ^ o exp (^\ % (aj - l)) |0). Though 
the coherent states are overcomplete, phase cancellations 
in Eq. (|27p leave the proper complete number basis at 
each t . 

By now-standard procedures [36, 38] we recognize that 
the fields <j> and <jr have somewhat different roles, with 
fluctuations in <f> about its mean value corresponding 
roughly to fluctuations in number, and those in <fr sam- 
pling moments of the generating functional \ip). Thus 
we expand the complex-conjugate coefficients <f>* of the 
(row) vector </>t at each time as <j)* = <pj + l at each t, with 

shorthand 4>\ = cfit + l, leaving <fi to be determined physi- 
cally. The resulting normalized generating functional has 
the path-integral representation 



(0|exp I J^oj ] |Vt) = / V^V^e'I^e^-^-^, 



in which the diffusive "Lagrangian" is 

L fa<t>) =4>- ^ + 0(0 + 1,0) . (29) 



The Liouville operator (IMf has induced a complex- 
valued function of fields (</>+ 1,<^>) through the sub- 
stitution a] — > <pj + 1, a,j — > 4>j. We will see that, up 
to care with signs and contours of integration that de- 
pend on what we wish to extract from this function, it 
behaves as the equivalent of a Hamiltonian for an equi- 
librium system, with a few structural differences charac- 
teristic of stochastic processes(elaborated also in 45]). 
The measure T>4>T>(j) in Eq. (|28|) is defined formally by 
the skeletonization procedure for insertion of the coher- 
ent states, but in practice is usually defined implicitly by 
perturbation theory in the diffusive Green's function [56j. 

Linearization of in either Eq. ([23|l or Eq. (|28|) (i.e. 
keeping only terms linear in <p) provides the natural 
expansion in independent collective fluctuations of the 
master equation and gives results for expectation values 
which are identical with the mean-field results presented 



earlier. In the field form (|29j) it further provides a conve- 
nient and intuitive background-field expansion, in which 
the backgrounds represent locally best-fit Poisson distri- 
butions with mean number equal to 4>*4>j for each com- 
ponent j. 



C. Structure of reaction-diffusion Lagrangians 

To make use of the form of f2, we introduce two projec- 
tion matrices onto the catalytic sites j = 0, J, P± with 
components {P+) J3 ' = ,j, and (P-)jj> = £j,o<5/,o> 

and linear-diffusion matrices 

1 



D 4 



-1 



1 

-1 



(30) 



D- = 



-1 
1 



-1 
1 



(31) 



corresponding to phosphorylation and dephosphorylation 
transitions, respectively. 

Noting that for 1 T the row vector of all ones, 1 T D± 
0. we can use dy or d> as it is convenient, to write 



1 



1 



1 



1 



1 + -tfP + 

q 



t£>+0. (32) 



We extract the overall iV-dependence of the action in 
the path integral by descaling time with the definition 
dz = dtq = dt/r, and rescaling the Lagrangian to a La- 
grangian density per particle, to write 



J dtL _ -N J dzL 



(33) 



with L = L/qN. If we similarly descale the field <fi ~^ 
<\> = <p/N, and the Liouvillian il — >• ft = Cl/qN, we have 
the Lagrangian density in terms of the natural coupling 
g = N/q: 



L = 6d z 



n 



(34) 



where 



Q (^,0) = (l + g^P-^j tfD_cj> 

l + g^P + 4^j^D+4>- (35) 
Note that the natural fields define the relative number 



operator 



rij/N = Uj, satisfying 
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Now not only are the fields 0' and expanded about dif- 
ferent backgrounds, comparable fluctuations of and 
correspond to fluctuations of 0^ and on scales differing 
by N, with large N defining the domain of perturbation 
theory. 

To expand the functional integral (|28p in Gaussian 
fluctuations, we further separate out mean values from 
the fields, introducing notation = + p (so putting 
0t = <p+ 1 T and (j>> = <ft + p) . Using a compact notation 
Cl'j for the tensor of i (^'-derivatives and j (^-derivatives 
of f2, the second-order Taylor expansion in p is exact: 



L = <j)-d z (j) + Vl 



ft»(*), 



(36) 



and Cl 2 is independent of 0t . 

The background = makes the first line of Eq. (|36|) 
vanish for general 0, and for more general we can 
expand in a classical background and perturbations, 
in which the linear order vanishes at that 0. The ip- 
linear term in the second line of Eq. (|36[) enforces a 6- 
functional if tp is rotated to an imaginary integration con- 
tour, and negative eigenvalues of Cl 2 (^cfij only soften the 

^-functional for their corresponding ip eigenvectors with a 
convergent Gaussian envelope. We handle these eigenval- 
ues in perturbation theory with a Hubbard-Stratonovich 
transformation [46| and a Langevin (auxiliary) field [38[ . 
We see below that in phases with no symmetry breaking, 

the eigenvalues of Cl 2 are all zero or negative. [57| 

Positive eigenvalues of Cl 2 (^j , of which one appears in 

the phase of symmetry breaking in this problem, require 
different treatment. They produce a divergent envelope 
for the ^-functional integral if tp is integrated along an 
imaginary contour, while a real contour for a p> eigen- 
vector does not enforce the expected (5-functional for the 
corresponding component of the diffusion equation. We 
expect, from experience with Euclidean field theories for 
reversible systems, that these eigenvectors signal the ex- 
istence of a continuous class of "approximate" stationary 
points generally termed instantons |35j . ip diverges ini- 
tially along a real contour, but for the appropriate joint 
background of 0, nonlinearities in the equations of mo- 
tion extend the divergence into a bounded trajectory of 
locally least J dzL, representing domain flips (a fluctua- 
tion that takes the system from one of the bistable phases 
to the other) in the symmetry-broken phase. The integra- 
tion over the unstable fluctuations of p are not handled 
in Gaussian perturbation theory about the static back- 
ground, but replaced (with a proper measure term) with 
the integral over all time-translates of the approximate 
stationary solutions. 



D. Symmetries and conservations 

Foregoing formal treatment of the convergence of 
Langevin perturbation theory and its regulation by ap- 
proximate stationary points [47J , we observe two impor- 
tant global symmetries of the theory which hold as field 
identities and also order-by-order in a large- N expansion. 
These are useful in numerically solving for approximate 
stationary points in low-dimensional examples. 

The classical equations of motion following from 
Eq. (|36l) and its equivalent expansion for = + 



d z 



fti (>,0) 



0. 



(37) 



(38) 



Both are 0(N°) , as L is defined in terms only of z, g, and 
descaled fields. The equivalent equations in terms of 0t 
and 0, resulting from shifts of the fields in the measure, 
generate Ward identities of the theory to all orders in N. 

The transformation $ — > e A 0^, — > e~ A at con- 
stant A is a symmetry of L at general (f>',<f>, whose as- 
sociated Nocthcr charge is number: d z (dfi ■ 0J = as a 

field equation. Time-translation is also a symmetry of L 

whose Noether charge is the potential: d z (Cl) = 0. Both 

of these follow immediately as properties of the classical 
solutions of Eq's. (|3T|) and (TJS|) . About backgrounds that 

are, or converge to, = 0, the constraints 0t • = 1 and 

Cl ^0~t, 0^ = specify a 2 J-dimensional subspace of field 

configurations in which all classical trajectories must lie. 
We further note that, due to the quartic form (f35|) . 



0-dJ = -0-A 1 (0"t,0) 

= -O(0~t,0) :Cl 2 (J). (39) 

The classical action over any stationary trajectory is then 

l(f~$) = ~\t :^ 2 (0) 

= -g (0P_0) (0£>-|) - g (0P+0) (0£ + 0\) . 

(40) 

The positive eigenvalues of Cl 2 (^j , which create di- 
vergent ip fluctuations if we include them in the expan- 
sion of Eq. (l36l) . correspond to trajectors that take L 
below the value (L = 0) of all classical (true) station- 
ary points. However, we will see that the nonclassical 
"approximate" stationary points of Eq's. ([3"T| and (|3"5)) 
produce strictly positive L, so that domain flips are sup- 
pressed relative to the persistence amplitude within do- 
mains in the symmetry-broken phase. This will be more 
transparent with the representation in terms of action- 
angle variables introduced in Sec. [V] 
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The background-field expansion to recover 
mean-field theory 



E. 



The relation (|28|) between the state and path integral 
representations of the master equation gives the expected 
first moment of the field 



) = (0| exp \J2 a j\ lV>t) = (n jt ) , 



(41) 



where the second angle bracket in Eq. (|4Tj) denotes expec- 
tation in the probability Pt(n). While not a field equa- 
tion (remember that <j)*j4>j is the combination extracted 
by fij), this relates the classical mean- field solutions to 
the stationary points of L. The classical solutions corre- 
spond to the subset of stationary points [H, H3] 



8L 
~d4> 







(42) 



= 



which solve Eq. (1371) at (ft = 1. These need not be time- 
independent, and include the full suite of classical diffu- 
sion trajectories. However, if the fij in the initial condi- 



tion ([26)1 are set to the steady-state values, they satisfy 
detailed balance under the ratio of catalytic rates corre- 
sponding to Eq. © (at A = in this case): 



x = = 



J 



1 



The fields themselves satisfy 

<f>j — x 3 (j>o, < j < J, 



(43) 



(44) 



per Eq. ([7]), and the remaining solutions for the coupling 
follow. 



IV. FLUCTUATIONS ABOUT STATIC MEAN 
FIELDS 

Fig. [S] shows the general character of timeseries for the 
population as represented by the number nj, in a phase 
with relatively strong symmetry breaking (g — 4.08 for 
J = 2. As a reference the phase transition occurs at 
g = 3). A timeseries is characterized by dense fluctua- 
tions about the mean value, in which nj remains near the 
mean-field value, punctuated by occasional large excur- 
sions that shift the mean. In this section we will consider 
the Gaussian-order approximation to the dense fluctua- 
tions about the mean. In Sec. [V] we return to qualita- 
tively different methods to handle the rare events which 
change mean population state. 



A. Diffusion eigenvalues and eigenvectors 

We compute the noise spectrum by further expanding 
Eq. (|3SJ) about = 0, defining <f) = 4> + 4 1 an d letting 
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FIG. 9: Simulation results for a population history, repre- 
sented by n.j (number of full-phosphorylated particles). We 
have used a system of 100 particles and a g-value of 4.08. 



4> be a constant solution to Eq. (|4"2")) . so that the linear 
term in (p vanishes. Using L ^0, (j^j — 0, the second-order 
expansion defines the Gaussian kernel, with the form 



L = (pDo(fi + LpDitp + h.o. 



(45) 



and higher-order terms (h.o.) are left for perturbative 
expansion. The diffusion kernel governing ij> in Eq. (|45[) 

is 



D a = d z + (l + g^ j D_ + (l + gj>j) D+ 

+ g (d_01 t P_ + D+4>1 T P+^ , (46) 
while the kernel controlling the constraint field tp is 



D 2 = | < | D 



V - /) ' <> /'j )+ transpose | 



(47) 

Remarkably, about general normalized solutions to 
Eq. 0) the kernel (l47l) has only two nonzero eigenvalues 
A± , with eigenvectors v± : 



D2V± = 2 A ± U ±- 



(48) 



We construct v± from convenient, orthonormal "center" 
and "edge" components, 



sinh (£) 



\j sinh [(J - 1)£] 
and zero otherwise, and 

(«e), 



' (^),l<i<J-l, 



>J=(J/2±J/2) V2cosh[(J-l)e] 
and zero otherwise. 



(49) 
(50) 
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A term that appears in the solution for the eigenvalues 
is abbreviated 



R(L) = \ li -Manh[(./ l.i^ianh (if). (.VI i 



in terms of which 



A± = 2 {±11 (£) - 1} cosh [(J - 1) f] 



' sinh(lg) 
sinh(^) 



(52) 



and the orthonormal eigenvectors 
«± = {^cV^ (0 ± 1 T (0 T l} ■ (53) 

For g < g c , only f = is consistent, and we get A + = 0, 

(54) 



J + l 



with eigenvector u_ = u e . This algebraic result em- 
phasizes the efficiency of expanding about Poisson back- 
grounds for weakly perturbed stochastic processes. The 
only deviation from Poisson which must be handled per- 
turbatively comes from a single mode of ip, whose fluctua- 
tions represent exchanges between no and rtj by Eq. (|50| . 
These are of course the noise in the catalytic rates that 
feeds back into the distribution as a whole. 



Q is (5-correlated in z with weight g/N, 



and drives the field 



(57) 



(58) 



through the inverse of Dq, acting 
on v-. In the symmetric phase A + = and this is all 
there is to the bare noise kernel; in the symmetry-broken 
phase we must still handle (by other means) the term 
ip ■ v+, which however remains orthogonal to the u_ in 
the Langevin term. Integration over ip in the symmetric 
phase produces 



4>z = \/-A- / dz'Go (z, z') 
Jo 



(59) 



as a field equation to Gaussian order, in which Go (z, z j 
is defined in terms of Eq. (|46|) by 



(60) 



Note that from Eq. [59] we see that \4> z ) — 0> which 
implies that there are no corrections to the mean-field 
result for the expectation value in the symmetric phase. 



D G (z,z') = S(z- z'). 



B. Hubbard-stratonovitch transformation about 
the symmetric phase 

Rather than complete the square in Eq. (|45t (a la On- 
sager and Machlup !48j), which is cumbersome for one 
eigenvector, we introduce into Eq. (|28|) an auxiliary-field 
representation of unity at each time: 



l=J\f V(e 



3 2g J s 



(55) 



in which is a time- and field-independent normaliza- 
tion. Shifting the auxiliary field £ (a symmetry of the 
measure), we introduce the physical Langevin field £ as 



C= C - flV-A- (ip ■ v-) 



(56) 



The net effect on L is the shift 
1 m 



2.9 



C + L 



C. Fluctuations about the symmetric order 
parameter 

As an example we compute to lowest order the fluctua- 
tions in the order parameter about the symmetric phase, 
where the diffusive Green's function is easy to compute 
in closed form. Application of the number operators 
first to the basis \n) and then to the coherent-states in 
Eq. (|28p yields the connected component of the variance 
(expressed in descaled fields) 



(("J - n Q ) 2 ^ - {n,j - n ) 2 

2N/JJ+T) = ' 



+N(J+l) 




From the field equation (|ST))) and the correlator (|5"5|) we 
obtain the <h noise 



V2 



-A_- 



' N 



dz'i 



[ -1 ••• 1 ] 

v 7 ! 



G [z,z') ■ 



(62) 
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and it remains only to compute the mode expansion of 
the diffusion kernel in D from Eq. (I46[) in the uniform 

background q> = 1/ (J + 1). 

The symmetric-phase Dq contains a symmetric 
linear diffusion matrix with endpoint corrections, 
so its eigenvectors v/- have components (ujt)j = 
Afk cos [9k + Kk (j ~ J/2)], with A/fc a normalization. The 
eigenvalues are immediate on the interior sites, 



In terms of these the modal expansion of the free 
Green's function is 



J+i 



Go (z,z') = Q(z-z')J2> 



fc=l 



(66) 



A, 



4 1 



.9 



J+l 



2 K fc 



and consistency of the interior with the endpoint correc- 
tions then determines n k and 9 k . 

Either sin ( J + 1) K k /2 = => K k = irk/ (J + 1), k 
even, and 6k = 0, or cos#fc = => = 7r/2 and 
solves the matching equation 



tan — 
1 



.9 



J + 1 + 2.9 



J+l 
tan Kfc. 



(63) 

Only odd-fc modes from Go couple to u_ in Eq. (I62[) . by 
symmetry, and the wavenumber sums from k > 3 are 
easily approximated with a two-dimensional A: integral. 
We note that for all these modes N 2 rj 2/ (J + 1), and 
the only values that contribute to the inner product come 
from j = +J/2, giving sin 2 (Kfc J/2) w 1. Thus the non- 
singular modes in the diffusion kernel provide a smooth 
background approximately linear in g. 



(64) 



Eq. (|6"4")l is regular for fc > 3 odd, and creates only a 
small wavenumber shift from the free diffusion solution, 
leaving n k w 7rfc/ (J + 1). The important mode for criti- 
cal behavior is k = 1 as 5 — > g c , where k± —> as 



1 



l) 2 V.9 



1 



1 

9c 



The leading contribution from the k = 1 mode oc- 
curs when it is present in both factors of Go- For small 
g — g c this mode is almost linear in j, with normalization 

Wf 2 -> K 2 Ej J £ 2 _.//2i 2 - Evaluating this singular term 
separately, with Eq. (1551 for the eigenvalue and Eq. 
for the limiting value of the wavenumber, and then com- 
bining with the background from the regular modes, we 
(65) obtain the approximation 



((nj-n ) 2 ) - (nj-Tio) 2 8 5 log (J+l) 



2N/ ( J + 1) 



7T (J+l) 



J+l 



M\9 



(67) 



It is convenient to separate the constant 

from the discrete sum for the k = 1 norm, because C — > 1 
at large J, but differs somewhat at the smaller J of more 
likely biological interest. 

The approximation (|6T[) to the closed-form mode ex- 
pansion for the variance of the order parameter is com- 
pared to numerical simulations in Fig. 1101 We have con- 
tinued the analytic expression through the critical point 
to show the peak, though the character of the modes 
rapidly changes as the distribution becomes skewed in 
the symmetry-broken phase. A similar mode expansion 
exists in this phase, but requires the relaxation eigenvec- 
tors and eigenvalues for asymmetric diffusion. We have 
not computed these in closed form, and do not pursue 
them numerically because in the symmetry-broken phase 
fluctuations quickly come to be dominated by the center- 



of-mass behavior we derive below in Eq. (|70p . 

Three main observations are important. First, the sin- 
gularity in the variance has the leading-order approxima- 
tion 

(("J - n of) - (nj-np) 2 const ~ 1 

2A7(J+1) (J + l) |S -Sc|' ( } 

comparable to that of the mean-field Ising ferromagnet, 
like the scaling of the order parameter. The variance 
has weight 1/ (J + l), because the lowest diffusive mode, 
corresponding to the average magnetization, is the only 
collective fluctuation participating in the phase transition 
near the critical point. 

Second, we see that the weak-coupling scaling of 

(nj — uq) 2 ^ — (uq — nj) 2 is that of Poisson noise for an 

average of NJ ( J + 1) particles per site. (This is shown 
in the inset of Fig. HUH 

Third - and the reason we do not pursue the low- 
order expansion for the symmetry-broken-phase two- 
point function - we see that for g 3> g c the variance 
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FIG. 10: Fluctuations in the order parameter scaled for g 
above and below critical. var(n') stands for the variance 
((no — n .j) 2 ') — (no — nj) 2 . Lines are leading-order expansion 
of Eq. (|67|l in fluctuations about the symmetric mean-field so- 
lution, continued through g c ; symbols from simulation, J + 1 
values and markers as in Fig. [5] Large panel shows conver- 
gence of var(n r ) to N/g at all J. Inset shows convergence of 
wxr(n') to Poisson result 2N/ (J + 1) as g — > 0. 



(inj — no) 2 ^ — (no — nj) 2 goes to a universal form N/g 

for any J. The independence of this scaling regime from 
J indicates that the particles interact with one end of 
the chain and the exogenous catalysis only, suggesting a 
strong-coupling limit. 

In addition, for large J, the center-of-mass equation ([3]) 
predicts (for I = P), in steady state for the symmetry- 
broken phase, that the variance will equal 



(n e )N 



1 

9 



[ n e) 

N 



(70) 



where n e signifies whichever of the end sites or J is 
occupied (which depends on which of the two bistable 
states is chosen). In Fig. [TT1 we plot the simulation data 
for the variance for J + 1 = 100 against the estimate 
from Eq. (1701 . For the value of (n e ), we simply take the 
value of the order parameter at the corresponding value 
of g. As we see, this explains the form of the fluctuation 
spectrum for large J very well. The fact that we are 
able to replace the order parameter \(nj — no) | by n e 
demonstrates that this scaling regime is independent of 
J as well and the particles only interact with one end of 
the chain all through the symmetry-broken phase. 



V. LARGE EXCURSIONS 

From the 1/g scaling of the Lagrangian term for the 
Langevin field in Eq. (|57[) and the presence of unsta- 
ble eigenvalues in the fluctuation kernel (p4"51 about static 
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FIG. 11: The fluctuation in the order parameter for J+l = 
100 and TV = 400 is plotted against the fit predicted by Eq. 

ran 



stationary backgrounds, we anticipate that perturbation 
theory containing fluctuations associated with domain 
flips will not converge [13], and that the rates for these 
large excursions will be computed as an essential singu- 
larity with respect to the perturbation expansion. The 
remarkable similarity to Hamiltonian dynamical systems 
created by this method for treating generating functions 
reduces the problem of estimating both escape trajec- 
tories and first passage times to that of identifiying the 
heteroclinic network [49|, [5(| in the associated dynamical 
system. 

We construct the expansion appropriate to the second- 
order transitions, beginning with the analysis of the ex- 
act stationary points corresponding to solutions of the 
classical diffusion equation from arbitrary initial condi- 
tions. From the structure of these solutions, we identify 
the "approximate stationary points" associated with do- 
main flips, and compute the trajectory and action for 
a low-dimensional example numerically. Comparisons to 
simulation suggest that this calculation correctly predicts 
the leading exponential dependence on particle number 
of the residence time in domains in the symmetry-broken 
phase. 



A. The expansion in semiclassical stationary points 

We define the stationary-point expansion of the path 
integral (1281) implicitly by the requirement that the resid- 
ual perturbation theory converge. In principle we must 
include not only the (generally unique) exact stationary 
point specified by the initial state \ipo), but also a suf- 
ficient set of "approximate" stationary points associated 
with states that converge exponentially fast toward |^o)) 
with respect to prediction of late-time observables. In 
this representation using generating functionals for prob- 
ability distributions over spaces of histories, a "stationary 
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point" refers to any full path ( (f>, $j satisfying the clas- 
sical condition that the linear variations of L (including 
time-derivative terms) vanish: 



dL „ dL n 
-~ =0; — = 0. 



(71) 



If <j> = 0, the solution of Eq. (|42[) solves the above to 
recover the mean- field solutions; in this section we relax 
the requirement <j) — to uncover a larger set of solutions 
which result in a nonzero value for L. 

Formally, then, the removal of the stationary-point 
contribute leads to the functional form 

(0|exp M*) = 

v<pv$ e - N ! H 1 - 1 ) e^-M , 

(72) 

where L denotes L((fc,(p\. As perturbative corrections 
scale as powers of 1/N, by Eq. (|62|) . to leading order 
we will treat L — L as a quadratic form in </? and <f> , of 



the form ((45]), now about more general time-dependent 
backgrounds. 

The formal sum ^ ~ ^ is properly a discrete sum in the 
number n G 0, . . . , oo of domain flips, of a time-ordered 
integral over their positions {z%, . . . , z n }. The integral 
is necessary (35[, because as the domain flips converge 
toward true stationary points, the fluctuation generated 
by time translation of any solution becomes a null eigen- 
vector of the functional determinant about that solution. 
(This is a variant on Goldstone's theorem (4(| , associated 
with the time-translation symmetry spontaneously hid- 
den by the instanton (45|.) This eigenvector is replaced 
by the integral (with a Jacobean), and the remaining 
functional determinant is a product of positive eigenval- 
ues, by construction. 

We will check that the transition times of the instan- 
tons are finite and that they converge exponentially to 
the static backgrounds, so that when they are improba- 
ble the dilute-gas sum is well defined. As we verify below, 
the classical solutions all have <f> = and zero action, and 

we denote by So the action N J dzL associated with a sin- 
gle instanton. Letting 1/Co denote the Jacobean relating 
the null eigenvalue to the measure for z-translation of the 
instanton, we recast the sum in Eq. (|72|l as 



(0|exp[]T J IV,) / ^,.-.,/ |/^3)|'e- , '/ i ( U )e« s -« (73) 



in which T denotes time-ordering in z of the positions of 
the instantons. The presence of n factors of 1/£q m the n- 
instanton determinant follows from the product structure 
of functional determinants and the wide separation of 
finite supports in z where the background differs from 
that of a steady state [35j|. 

Like the computation of the energy shift in the equi- 
librium double-well problem, we see that observables re- 
lating to persistence within a domain will receive con- 
tributions from the even terms in the sum over n, while 
those relating to domain flips will receive contributions 
from the odd terms. The likelihood of persistence decays 
exponentially in z at early times with rate 

rflip = ^e~ S °. (74) 
Co 

The computation of the instanton action So, which is 
responsible for the leading exponential dependence of rm p 
is most easily carried out within the complete analysis of 
the semiclassical stationary points, beginning with the 
classical diffusion solutions. 



B. Action-angle variables and the structure of the 
Hamiltonian 

The equations of motion Q37I38I) in (f> and <f> do not 
directly give the evolution of the physical particle num- 
bers, or efficiently use the symmetries of Sec. IIII D[ To 
do both, it is convenient to transform the background 
fields as . = e (7j , and <j>j = Vje~ aj . Uj is then the semi- 
classical approximation to the relative number (n.j) /N. 
This change of variables is equivalent to an action- angle 
transformation in classical mechanics [51| , and we check 
in Ref. [45[ that as well as producing a more convenient 
form for the action, it leads to the correct measure for 
fluctuations. The Lagrangian (|34p retains a simple ki- 
netic term, up to a total derivative: 

L = a ■ d z v + £l(a,v) , (75) 

and the equations of motion in the new variables become, 
respectively, 
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and 



dVj 



(77) 



With the interpretation of the number field v as a po- 
sition, and er its canonically conjugate momentum, — fl 
becomes the correctly signed Hamiltonian for classical 
solutions. The conservation law dCl/dz — is math- 
ematically a conservation of energy, but the particular 
value O = associated with all stationary points initi- 
ated by classical distributions is a distinctive feature of 
this stochastic-process application of Hilbert-space meth- 
ods. 

The global symmetry whose Noether charge is total 
number becomes immediate in action-angle variables. 
Defining 



+ } (crj - a) d z Vj + tl (a, v) 



3=0 



(79) 



in which a multiples the z-derivative of conserved total 
number in the first line, and only differences Ci — aj ap- 
pear in either the kinetic term or the f2 of the second 
line. 

To expose the structure of the associated dynamical 
system, and to make the terms in it readily visualizable 
from the mean-field diffusive solutions, we perform a final 
transformation by introducing the log ratio of particle 
fluxes between sites j and j + X, 



1 V 



(78) 



3=0 



the Lagrangian becomes 



L 




r 3 + l-3 = 2 |U " 



Q + gvj)vj 



(1 + gv ) Vj+i 



(80) 



In terms of these the er-dependence of Cl may be simplified 
to read 



j-i 

Cl (cr, z/) = 2y/ (1 + gvp) (1 + gvj) {^Jv j+ iVj [cosh (r j+1 j) - cosh (aj - a J+x - r j+lij )]} . (81) 

3=0 

The one-dimensional geometry we have assumed for the graph of phosphorylation and dephosphorylation transitions 
makes possible the definition of a function a(v) at each number configuration, satisfying 



a j+1 (v) - a,j(v) = -r j+1J , 



(82) 



(up to a convention for specifying a at each v, which we may take to be arbitrary). a(v) is a reference point for the 
canonical momentum coordinate cr, and a — a(v) functions as the kinematic momentum for this system. 
We may see this by introducing the Hamiltonian potential function 



J-i 



V(v) = 2yJ (1 + gv ) (1 + gvj) ^ v^+i^ [cosh (r j+1J ) - 1] 

3=0 



in terms of which 



J-i 



(cr, v) = 2^/(1 + gvo) (1 + gvj) ^ V^j+^j { cosh ( a ~ a (v))j ~ ( a ~ a i v ))j+l -l}-V(v) 



(83) 



(84) 



3=0 



To leading order the explicit cosh in Eq. (1541 is sim- 
ply a quadratic form in a — a(v), with a matrix of in- 
verse masses defined by the remaining square-root terms. 
When a — a(v) = 0, Eq. (jTo) shows that d z v = 0, verify- 
ing the interpretation of a — a(v) as the kinematic mo- 
mentum. Furthermore, if this kinematic momentum van- 



ishes at any local minimum of V(v), Eq. (|77|) shows that 
d z a = 0, hence d z (a — a(v)) = 0. The local minima sat- 
isfy V(v) = and are attained only when each rj+ij = 
independently, because the cosh terms in Eq. ([531 are 
never negative. These are of course exactly the (stable 
and saddle-point) mean-field solutions with particle ex- 
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change between adjacent sites obeying detailed balance. 
We identify them in graphs below as the fixed points of 
the classical diffusion equation. 

It can be shown [4o| that, as long as the quadratic ex- 
pansion in a is a good approximation, and as long as the 
effective mass terms implicit in Eq. (|84[) are not a strong 
function of v (which we will verify), all stationary points 
of the action closely approximate ordinary mechanical 
trajectories in the potential —V(v), with position coordi- 
nate v and kinematic momentum a — a(y). For the classi- 
cal solutions (7 = 0, shown in Fig. [T^l the unbounded tra- 
jectories are those that originate in non-equilibrium ini- 
tial conditions and converge exponentially slowly on the 
saddle or stable fixed points. The two bounded trajecto- 
ries, between the saddle and either stable fixed point, 
travel along the saddle path of the potential, —V(v), 
which is bounded above by and unbounded below, and 
make up part of the heteroclinic network (49l . |50] ] of the 
associated hyperbolic system. 




horizontal displacement on the unit simplex 



FIG. 12: Classical flowfield of the particle numbers Uj for 
J — 2, shown on the simplex Yl'j-o v i ~ ^' Upper left corner 
is no = 1, upper right corner is ri2 = 1, and bottom is ni = 1. 
Dots are a selection of initial conditions, and classical diffusion 
solutions are the flowlines emanating from them. The uniform 
distribution Uj = 1/3 (open circle) is the saddle point, and 
the stable fixed points (heavy dots) are attractors. 



For ordinary mechanical Sow, we know that the full 
heteroclinic network consists of trajectories running both 
ways between the stable and saddle fixed points. If this 
were a purely mechanical system, the reverse trajecto- 
ries would be strict time-reversal images of the bounded 
classical diffusion trajectories. (Note that an exact rever- 
sal (<7j — <7j+i) — > 2rj + ij — (<7j — o'j+i) would also leave 
Cl = 0.) Here, a small ^-dependence of the effective mass 
terms causes them to differ slightly from each other and 
from the saddle path over —V(v). 



In Fig. [T3l we directly compute the trajectory of the 
reverse bounded path, by integration along the saddle 
instability of the equations of motion (|76|77|) . The fact 
that it nearly retraces the classical diffusive direction of 
slowest flow checks the approximation that both trajec- 
tories are dominated by the potential —V{v) itself. 




horizontal displacement on the unit simplex 

FIG. 13: Flowfield of the instanton solution in J = 2, for 
values £ = 1, so g m 4.0862. Fine lines are classical diffusion 
solution from Fig. 1121 and bold line is the instanton. 



For an instanton in a classical equilibrium field theory, 
the conjugate and the kinematic momentum would be 
the same quantity. Both forward and reverse trajectories 
along the saddle path in the potential —V{v) would have 
locally minimum but non-zero action, and in that sense 
both would be "non-classical" trajectories (35|. The dis- 
tinctive feature of the path integrals associated with mas- 
ter equations of the form we have considered here is the 
offset cr(y) from the canonical momentum that appears 
in the kinetic term of the Lagrangian (|T5j) to the kine- 
matic momentum. This offset is responsible for S = for 
all diffusion solutions, including the bounded trajectories 
from saddle to stable fixed points, and it approximately 
doubles the value of the Lagrangian (fTSj) along the re- 
verse trajectories. The value of this Lagrangian, along 
the numerically determined path of Fig. [T31 is shown in 
Fig. HH Like the Lagrangian for a classical problem, 
it is positive definite, and approximates the WKB inte- 
gral for barrier escape, except with an extra factor of 2: 
S w 2N j Vdv T mdisy/2V(is), where dv is a length ele- 
ment on the coordinate v, and m is the matrix of effective 
mass values implied by Eq. (|84[) . This approximate form 
follows simply from the nearly time-reverse character of 
escapes versus classical paths of slowest-diffusion, and 
may be derived from the original Gaussian-order approx- 
imation to such escapes by Onsager and Machlup [48j . 
Readers seeking a systematic derivation, including the 
Onsager-Machlup small-fluctuation approximation, may 
find these in Ref. [H| or Ref. [Uj]. 
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FIG. 14: Linear-scale plot of the Lagrangian density per 
particle L of the instanton trajectory. A logarithmic plot 
of the same quantity (not shown) demonstrates exponential 
decay toward zero at early and late times. 



For the parameter values of Fig. [T3J the integral J dzL 
under the curve of Fig. [TJ] converges to a value near 
0.0206. Fig. [15] compares numerical estimates of the res- 
idence time in this model, inverse to the rate rfi lp of 
Eq. (|74[) . to particle number N. The slope of the log- 
arithm of 1/Vflip should be dSo/dN = J dzL, up to cor- 
rections decaying as 1/iV, and we observe quantitative 
agreement with the numerical estimate of the instanton 
to - 10%. 

Although we do not pursue the analytic forms in this 
paper, the dependence of the coefficient of N (giving de- 
cay times [Hj]) on number of sites J and the distance from 
criticality g — gc may also be found from simulations. 
Fig. I16l svmmarizes these numerical results, showing that 
the dependence on small-integer J is roughly linear, and 
the dependence on g — gc < 1 is weakly nonlinear. 

C. Summary of semiclassical results 

We have seen that the classical "action" for this 
reaction-diffusion theory, once constructed, yields quite 
nicely the two extremes of behavior of interest in coop- 
erative intermolecular phase transitions. The classical 
stationary points coincide exactly with the usual mass- 
action differential equations. Corrections to these from 
cubic and higher-order fluctuation effects are readily in- 
corporated (for an example, see Ref. [Hj]), but to the 
resolution of our simulations, we cannot identify a need 
for such corrections at these parameter values, so we have 
not pursued them. 

The nonclassical stationary trajectories are the projec- 
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FIG. 15: Plot comparing residence times to particle number 
for the parameters of Fig. 1141 Slope of numerical estimates, 
which should correspond to dSo /dN — J dzL is best fit by 
0.023, while direct integration under the curve of Fig.rTilvields 
0.0206. 



tion from this (2 J + 2)-dimensional configuration space, 
onto the one-dimensional path most likely to destabilize 
the symmetry-broken phase, which closely approximates 
the path of slowest diffusive correction. The methods 
shown here therefore provide a compact and convenient 
way to estimate the escape trajectories and first-passage 
times for even quite richly structured nonlinear diffusion 
processes of this kind. These methods, originally de- 
veloped for applications to reaction-diffusion theory, are 
increasingly finding applications in epigenetics [ijg, IH| 
and systems biology [401 ] where particle numbers may be 
small, making fluctuation effects important, while at the 
same time the structure of the state space remains com- 
plex to describe. 

The reduction to a onc-dimcnsional system was as- 
sumed given, in the treatment in Ref. [l3j of a switching 
system comparable to ours; we have shown here a system- 
atic approach to estimating such escape trajectories. We 
have also verified that the action of the instanton, eas- 
ily numerically integrated once the trajectory is found, 
produces both the correct N scaling, and good quanti- 
tative agreement, with the domain residence times. We 
expect that, while the treatment of the functional deter- 
minant will be more difficult for metastable domains in 
first-order transitions, the classical-level analysis compa- 
rable to ours will be similar, and roughly as effective. 

VI. DISCUSSION AND CONCLUSIONS 

We have tried to idealize in a reasonable way a large 
class of biomolecular signal transduction systems, and 
to apply the most complete formalism available to de- 
compose and quantitatively estimate their properties as 
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FIG. 16: Plot of expected residence times r from simulations 
in the symmetry-broken state as functions of J, N, and dis- 
tance of the coupling strength from its critical value, g — gc- 
(a): logr vs. N for J £ {2,..., 5}. In simulations I = P, 
g = 7V/V IP is held constant, and g ~ gc is held equal to 
1. Green curve for J = 2 and g = 4.0862 corresponds to 
the value £ = 1, which we have computed analytically, corre- 
sponding to the coefficient dlogr /dN shown in Fig. 1151 (b): 
logr vs. g — gc for values J £ {2, . . . , 9} (from bottom to 
top) at TV = 50. 



switches. Our results thus combine a number of techni- 
cal advances in recognized domains, with several concep- 
tual insights relevant to the robustness and evolution of 
devices. Of necessity in a short treatment, our idealiza- 
tions of feedback and single timescales for all microscopic 
motion have abstracted away from some important prob- 
lems of connection with phenomenological models of real 
regulatory protein systems. 



A. Technical advances 

The operator treatments of gene-expression switch- 
ing Q extended traditional mass-action models to in- 
clude perturbative noise from first principles. We have 
further extended the operator methods to a path-integral 
treatment, which adds an intuitive and computationally 
tractable approach to large deviations. 

We have demonstrated that expansion of weakly non- 
linear stochastic processes about Poisson backgrounds 
leads to very efficient perturbative schemes for correct- 
ing the full probability distribution (not very surpris- 
ing in retrospect), and that in our particular idealized 
model, the entire noise spectrum is driven by a single 
bare Langevin field (perhaps somewhat more surprising, 
and not noticed before). 

Finally, we have shown that the nonlinear projection 
of the full master equation onto the dominant trajec- 
tory participating in domain flips approximately, but not 
exactly, reverses the unique trajectory of slowest diffu- 
sive correction in the classical flow. We have recovered 
the exponential in N characteristic of extensive large- 
deviations scaling [43| , and shown how to estimate the ex- 
act coefficient to refine the bounds of order unity that are 
conventionally (and usually correctly) assumed in pure 
scaling arguments [l3j . 



B. Biological insights 

The most concrete of our results for biologists seek- 
ing to understand the function of signal-transduction cas- 
cades and switches is that particle number (N), as well as 
exogenous kinase (I) and phosphatase (P) numbers, can 
be used to control the onset of switching, and in cases 
of asymmetric topology, also the preferred domain of the 
switch. The control through N offers a feedback from 
gene expression into the function of the cascade, which 
apparently has not been considered earlier. 

We have demonstrated through the mode expansion 
in the neighborhood of the phase transition, that only 
the lowest-eigenvalue collective fluctuation of the diffu- 
sion operator induces the instability to symmetry break- 
ing, and scales the divergence in the noise spectrum. For 
large J, we can also predict fluctuations in the symmetry- 
broken phase because of a simplification induced by the 
fact that the system senses only one boundary in the en- 
tire symmetry-broken phase. This enables us to simplify 
an exact equation for the center-of-mass of the system to 
predict fluctuations. 

More abstractly, we have distinguished a mechanism 
for switching based purely on population-level polariza- 
tion of the protein pool, from mechanisms which depend 
on limiting one or more transition rates through restric- 
tions on catalytic kinetics. Polarization-based mecha- 
nisms make the function of the switch dependent on its 
catalytic topology and concentrations, but not on kinetic 
factors, a separation that has been proposed as a route to 
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modularity. We hope that such distinctions can at some VII. ACKNOWLEDGEMENTS 

point be incorporated in evolutionary models that make 
quantitative use of the structure of phcnotypically neu- 
tral networks, where we hope they will explain at least 
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